\documentclass[11pt]{article}
\usepackage{amsmath,amsfonts,afterpage}
%\usepackage{showlabels}
\usepackage[pdftex]{graphicx,color}
\usepackage{listings}
\newcommand{\normal}[2]{{\cal N}(#1,#2)}
\newcommand{\normalexp}[3]{ -\frac{1}{2}
      (#1 - #2)^T #3^{-1} (#1 - #2) }
\newcommand{\La}{{\cal L}}
\newcommand{\fnom}{\tilde f}
\newcommand{\fhat}{\hat f}
\newcommand{\COST}{\cal C}
\newcommand{\LL}{{\cal L}}
\newcommand{\Prob}{\text{Prob}}
\newcommand{\field}[1]{\mathbb{#1}}
\newcommand\REAL{\field{R}}
\newcommand\Z{\field{Z}}
\newcommand{\partialfixed}[3]{\left. \frac{\partial #1}{\partial
      #2}\right|_#3}
\newcommand{\partiald}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\argmin}{\operatorname*{argmin}}
\newcommand{\argmax}{\operatorname*{argmax}}
\newcommand\norm[1]{\left|#1\right|}
\newcommand\bv{\mathbf{v}}
\newcommand\bt{\mathbf{t}}
\newcommand\Vfunc{\mathbb{V}}
\newcommand\Vt{\mathbf{V}}
\newcommand\vexp{V_{\rm{exp}}}
\newcommand\texp{T_{\rm{exp}}}
\newcommand\cf{c_f}
\newcommand\cv{c_v}
\newcommand\fbasis{b_f}
\newcommand\vbasis{b_v}
\newcommand\tsim{{\mathbf t}_{\rm{sim}}}
\newcommand\DVDf{\partiald{\Vt}{f}}
\newcommand\Lbb{\mathbb{L}}
\title{Notes on EOS Estimation Code}

\author{Andrew M.\ Fraser}
\begin{document}
\maketitle
\begin{abstract}
  While source for these notes reside at
  https://code.google.com/p/metfie/like\_gun/notes.tex and the
  illustrations are in the same directory, I use the same ideas for
  work on the 9501 model on moonlight.
\end{abstract}

\section{Notation for and Implementation of Basic Functions}
\label{sec:basic}

Given an experimental sequence of measured pairs of velocity and time
values, $(\vexp, \texp)$ and an initial model, $\fnom$, I describe how
to estimate a new model, $\fhat$.

\subsection{Notation}
\label{sec:basic_notation}

\begin{description}
\item[$(\vexp, \texp)$:] An experimental sequence of measured pairs of
  velocity and time
\item[$\bv$:] A sequence of model velocities
\item[$\Vt$:] A map from times to model velocities, eg, $\bv =
  \Vt(\texp)$
\item[$\tsim$:] A sequence of closely spaced sample times at which to record
  simulated position and velocity for constructing $\Vt$.
\item[$f$:] An EOS function
\item[$\Vfunc$:] An expensive procedure that maps an EOS function to a
  function that maps time to velocity with $\Vt = \Vfunc(f)$
\item[$\cv,\vbasis$:] Vectors of spline coefficients and basis functions
  that define a $\Vt$
\item[$\cf,\fbasis$:] Vector of spline coefficients and basis functions
  that define an EOS
\item[$D$] Matrix with elements $D[j,k] = \partiald{\cv[j]}{\cf[k]}$
\item[$B$] Matrix with elements $B[i,j] = \vbasis[j](\texp[i])$, ie,
  the exactly linear map from $c_v$ to simulated data.
\item[$\DVDf$:] The derivative of $\Vt$ with respect to $f$; expensive
  to evaluate.  $\DVDf (\texp)$ is represented by the matrix $(BD)^T$.
\item[$\epsilon$] The vector of differences between the measured
  velocities and the simulated velocities at the times of the
  measurements.
\end{description}

\subsection{Implementation}
\label{sec:basic_implementation}

\begin{description}
\item[$\Vfunc(f)$:] Run a simulation and record a sequence of
  velocities $\bv$ at times $\tsim$.  Then fit a
  spline to $(\bv, \bt)$ to obtain $\Vt\iff \{cv,\vbasis\}$.
\item[$\Vt$:] Call the spline evaluation method to get $\bv = \Vt(\texp)$
\item[$f$:] A spline with coefficients $\cf$ and basis functions
  $\fbasis$
\item[$\DVDf$:] Evaluate $\Vfunc(f)$ for $N(\cf)+1$ different sets of
  coefficients $\cf$, and use finite differences to get a matrix with
  elements $D[j,k] = \partiald{\cv[j]}{\cf[k]}$.  Then
  \begin{equation*}
    \partiald{\Vt(T)}{\cf[k]} = \sum_j D[j,k] \vbasis[j](T).
  \end{equation*}
\end{description}

\section{Priors, Likelihood and Optimization}
\label{sec:opt}

\subsection{Minimum Squared Error}
\label{sec:minsq}

Before calculating the step for the full log a posteriori probability,
I calculate a least squares step (maximizing the likelihood if all
$\sigma_v[i]$ are equal to each other).  Thus, I wish to find $d$ that
minimizes
\begin{equation}
  \label{eq:ssq}
  \tilde S(d) = \sum_{i=1}^{N_{\rm exp}} (\vexp[i] - \Vt_{f+\delta}(\texp[i]))^2,
\end{equation}
where $f+\delta = \sum_k (\cf[k] + d[k])\fbasis[k]$, and 
\begin{equation}
  \label{eq:taylor}
  \Vt_{f+\delta}(\texp[i]) = \Vt_f(\texp[i]) +
  \sum_j \sum_k \partiald{\cv[j]}{\cf[k]}d[k]\vbasis[j](\texp[i])
  +\text{ HOT}.
\end{equation}
Defining
\begin{equation}
  \label{eq:epsilon}
  \epsilon[i] \equiv \vexp[i] - \Vt_f(\texp[i])
\end{equation}
and dropping HOT, I write
\begin{align}
  S(d) &\equiv \tilde S(d) - \text{HOT} = \sum_i \left( \epsilon[i] -
  \sum_j \sum_k \partiald{\cv[j]}{\cf[k]}d[k]\vbasis[j](\texp[i])
  \right)^2 \nonumber \\
  \label{eq:sd}
  &= (\epsilon - BDd)^T(\epsilon - BDd)
\end{align}
and
\begin{align*}
  \partiald{S}{d[l]} &= -2 \sum_i \left( \epsilon[i] -
  \sum_j \sum_k \partiald{\cv[j]}{\cf[k]}d[k]\vbasis[j](\texp[i])
  \right)
  \sum_j \partiald{\cv[j]}{\cf[l]}\vbasis[j](\texp[i])\\
  \partiald{S}{d} &= -2\left(BD\right)^T\left(\epsilon - BDd\right).
\end{align*}
Thus I seek $\hat d$ that solves
\begin{equation}
  \label{eq:dhat}
  \left(BD\right)^T\epsilon = \left(BD\right)^T BD \hat d.
\end{equation}

\subsection{A Posteriori Probability}
\label{sec:app}

I use Gaussians with constant diagonal covariances for priors and
likelihood as follows:
\begin{align}
\vexp,\texp | f+\delta &\sim
\normal{\vexp-\Vt_{f+\delta}(\texp)}{\Sigma_v} \\
\cf+d &\sim \normal{\cf+d - c_{\fnom}}{\Sigma_f}
\end{align}
Up to irrelevant constants $K$ the log probabilities are:
\begin{align}
  L_c(\cf+d) & \equiv \log(K_2) + \log \left(\Prob(\cf+d) \right)\\
  & = \normalexp{\cf+d}{c_{\fnom}}{\Sigma_f}\\
  L_v(\Vt_{f+\delta}) & \equiv \log(K_1) + \log
  \left(\Prob(\vexp,\texp | f+\delta) \right)
  \nonumber \\
  & = \normalexp{\vexp}{\Vt_{f+\delta}(\texp)}{\Sigma_v} \\
  &= \normalexp{\epsilon}{BDd}{\Sigma_v}.
\end{align}

Given $f$ (in terms of $\cf$) and the functions $\Vt$, $\DVDf$, etc.\
that depend on it, an optimization step consists of ignoring the
\emph{higher order terms} (HOT) and solving for the vector $d$ that
maximizes
\begin{equation}
  \label{eq:L}
  \Lbb(\cf,d) = L_c(\cf + d) + L_v(\Vt_{f+\delta}(\texp)).
\end{equation}
Differentiating, I find
\begin{align*}
  \partiald{\Lbb(\cf,d)}{d} &= (c_{\fnom} - \cf - d)^T
  \Sigma^{-1}_f + (\epsilon - BDd)^T \Sigma_v^{-1}BD.
\end{align*}
Thus I seek $\hat d$ that solves
\begin{equation}
  \label{eq:dmap}
  (BD)^T\Sigma_v^{-1}\epsilon + \Sigma_f^{-1} (c_{\tilde f} - \cf) = 
  \left((BD)^T\Sigma_v^{-1}BD + \Sigma_f^{-1} \right) \hat d.
\end{equation}

\subsection{The Code \emph{calc.py}}
\label{sec:code}

In the following methods of class \emph{GUN} which implement the
optimization, the values of the first and last 5 spline coefficients
of $f$ are not varied.
\begin{description}
\item[set\_t2v] Run a simulation to build a spline for mapping times to
  velocities and save as self.$t2v$.  The result for the initial EOS
  appears in Fig.~\ref{fig:vt_test}.
  \begin{figure}
    \centering
    \resizebox{\columnwidth}{!}{\includegraphics{vt_test.pdf}}  
    \caption{Velocity as functions of time.  Traces for the
      \emph{experimental} data, a simulation with the initial EOS and
      the errro.}
    \label{fig:vt_test}
  \end{figure}
\item[set\_D] Use finite differences to calculate dv/df in terms of
  spline coefficients and save as self.$D$ with
  $D[j,i] = \frac{\Delta c_v[j]}{\Delta c_f[i]}$.  See
  Fig.~\ref{fig:D_test} and note that $i$ and $j$ range over variable
  components of $c_v$ and $c_f$ (excluding the last 4 which are always
  0).
  \begin{figure}
    \centering
    \resizebox{\columnwidth}{!}{\includegraphics{D_test.pdf}}  
    \caption{The matrix $D$.  Each trace represents
      $\frac{\partial c_v[i]}{\partial c_f[j]}$ for a different $i$.}
    \label{fig:D_test}
  \end{figure}
\item[set\_Be] Map experimental velocities v and times t to the
  following:
  \begin{align*}        
    \text{self.}ep[i] &= v[i] - t2v(t[i]) \text{ Difference between simulation
                and data} \\
    \text{self.}B[i,j] &= b_j(t[i])     
  \end{align*}
  Notes: $b_j$ is the jth basis function for the t2v spline.  $B$ is
  the derivative of the simulated velocity function at the
  experimental sample times wrt the coefficients of the spline fit
  that implements that function, and $BD_{i,j}$ is the derivative of
  $v(t[i])$ wrt $c_f[j]$.  See Fig.~\ref{fig:BD_test}.
  \begin{figure}
    \centering
    \resizebox{\columnwidth}{!}{\includegraphics{BD_test.pdf}}  
    \caption{The matrix $B\cdot D$.  Each trace represents
      $\frac{\partial v(t)}{\partial c_f[j]}$ for a different $j$.}
    \label{fig:BD_test}
  \end{figure}
\item[func] From \eqref{eq:sd} the objective function is:
  \begin{equation*}
    S(d) = \left| \epsilon - BD\cdot d \right|^2.
  \end{equation*}
  The difference between the simulated and measured velocities is
  $\epsilon$, and $d$ is the change in the vector of spline
  coefficients.  $BD$ is the derivative of the simulated velocities
  wrt to $d$.  Note that $d\in \REAL^{n_c}$ where $n_c = len(c_f)-4$.
\item[d\_func] Derivative of the objective function,
  \begin{equation*}
    \partiald{S}{d} = -2(BD)^T(\epsilon - BD\cdot d)
  \end{equation*}
  See Fig.~\ref{fig:d_func_test}
  \begin{figure}
    \centering
    \resizebox{\columnwidth}{!}{\includegraphics{d_func_test.pdf}}  
    \caption{The derivative of the objective function $S$ (see
      Eqn.~\eqref{eq:sd}) with respect to $c_f[j]$ plotted against $j$.
      The first trace $dS_0$ is for the nominal EOS, and $dS_1$ is for
    the EOS that solves \eqref{eq:sd}.  The corresponding EOS
    functions appear in the lower plot.}
    \label{fig:d_func_test}
  \end{figure}
\item[free\_opt] Solve Eqn.~\eqref{eq:sd} for $d$ to obtain $\hat d$
  using the SVD solver numpy.linalg.lstsq.  The solution appears in
  Fig.~\ref{fig:d_hat_test}.
  \begin{figure}
    \centering
    \resizebox{\columnwidth}{!}{\includegraphics{d_hat_test.pdf}}  
    \caption{The function $\hat d$ that solves Eqn.~\eqref{eq:sd} for $d$.}
    \label{fig:d_hat_test}
  \end{figure}
\item[constraint] Return the vector of constraint function values:
  \begin{description}
  \item[$r_{i} = \left. \frac{d^2 f}{d t^2} \right|_{t_i}~\forall i
    < n_t$:]
    The second derivative of $f$ wrt to $t$ at all $n_t$ unique knot
    values for $t$.
  \item[$r_{n_t} = - \left. \frac{d f}{d t} \right|_{t_{n_{t}-1}}$:]
    Minus the derivative of $f$ at the last unique knot.
  \item[$r_{n_t+1} = f(t_{n_{t}-1})$:] The value of $f$ at the last
    unique knot.
  \end{description}
  If these constraints are satisfied, ie, $r_i>=0~\forall i<n_t+2$,
  then the function $f$ is positive, monotonic and convex over its
  entire domain.  
\item[calc\_d\_constraint] Calculate the derivative (wrt to c) of the
  constraint function, store result as self.d\_constraint and return
  it.
  \begin{description}
  \item[$r_{i,j} = \frac{d}{d c_j} \left. \frac{d^2 f}{d t^2}
    \right|_{t_i}~\forall i < n_t$:]
    The derivative wrt $c_f$ of the second derivative of $f$ wrt to
    $t$ at all $n_t$ unique knot values for $t$.
\item[$r_{n_t,j} = - \frac{d}{d c_j} \left. \frac{d f}{d t}
    \right|_{t_{n_{t}-1}}$:]
  \item[$r_{n_t+1,j} = \frac{d}{d c_j} f(t_{n_{t}-1})$:]
  \end{description}
\item[opt] Do a constrained optimization step.
  The goal of \emph{opt} is to find the change $d$ in spline
  coefficients $c_f$ that yield a function $f(c_f+d)$ that minimizes $F(f)$.
  The key to \emph{opt} is the following call to
  scipy.optimize.fmin\_slsqp
  \lstinputlisting[language=Python, firstline=366, lastline=376]{calc.py}
\end{description}

\subsection{Numerical Results}
\label{sec:numerical-results}

\newcommand{\freq}{k} %
Figures \ref{fig:d} and \ref{fig:fve} illustrate iterative maximum
likelihood fitting of an EOS.  The procedure starts with the following
initial function
\begin{equation}
  \label{eq:initial}
  \fnom(x) = \frac{C}{x^3},
\end{equation}
and the iterations approach the \emph{actual} function
\begin{align}
  \label{eq:actual}
  f(x) &= \fnom(x) + \frac{2 e^{-\frac{x^2}{2w^2}}
         \fnom(x_0)}{\freq}  \sin(\freq x)\\
   \text{ where }C &= 2.56\times 10^{10} \nonumber \\
  \freq &= 0.2 \nonumber \\
  w &= 0.2, \nonumber
\end{align}
which was used to generate the simulated data.
\begin{figure}
  \centering
    \resizebox{\columnwidth}{!}{\includegraphics{fig_d.pdf}}  
    \caption{Illustration of the derivative of the velocity function
      with respect to the pressure function, $\frac{d v}{d p}$.  The
      second and third columns correspond variations of the
      coefficients of basis functions that differ by a factor of ten.
      The lower left plot indicates that nonlinear part of the
      response is a factor of a thousand smaller than the linear
      response.  }
  \label{fig:d}
\end{figure}
\begin{figure}
  \centering
    \resizebox{0.8\columnwidth}{!}{\includegraphics{fig_fve.pdf}}  
    \caption{Maximum likelihood estimation of $f$.  The \emph{true}
      EOS appears as \emph{actual} in the upper plot, and the
      optimization starts with the \emph{initial} and ends with
      \emph{fit}.  The corresponding velocity as a function of
      position appears in the middle plot, and the sequence of errors
      in the velocity time series for each step in the optimization
      appears in the lower plot. }
  \label{fig:fve}
\end{figure}

\section{Comparing Model Classes}
\label{sec:classes}

We want to test the hypothesis that splines are better for building
models of an EOS than polynomials.  We hope to compare the bias and
variance of two sequences of estimators, one that uses splines and one
that uses polynomials.  We will study estimators for quantities of
interest that are functionals of the EOS that differ from the
functionals that correspond to the measurements.  In previous work, we
have analyzed the functionals for muzzle time and muzzle velocity for
an ideal gun, and in the sections above we have described an analysis
of velocity measurements that depend on time.  We hope that those
quantities of interest and those measurements are sufficiently
different in that the measurements are given as functions of time and
the quantities of interest are required as functions of position.

Here is the agenda:
\begin{enumerate}
\item Design two classes of EOS models, ie, pressure as a function of
  volume along an isentrope.  One class will be splines with $n$ knots
  equally spaced on a log scale, and the other will be $n^{\text{th}}$
  order polynomials.
\item Write code that draws realizations of experimental data based on
  an EOS that is outside of both model classes for any finite $n$.
\item Write code that estimates the coefficients of polynomial models
  of the EOS.
\item Calculate (or estimate with simulations) sequences (indexed by
  $n$) of bias and variance pairs of estimators for quantities of
  interest using the two classes of models.  We hope to find that the
  curve defined by the pairs for splines lies below the curve defined
  by pairs for polynomials.
\end{enumerate}

\subsection{Alternative Measurements}
\label{sec:alternative}

If we need measurements that are more different from the quantities of
interest, we have could use the following simplification of
measurements of cylinders.  Suppose that all of the material in a long
cylinder detonates in an instant and that the density after that is
uniform inside.  Thus the function $\left\{r(t): t>0 \right\}$ is a
sufficient specification of an experimental result.  Letting
$r(0) = r_0$ denote the initial radius the work done by the gas in
expanding to $r$ is
\begin{equation*}
  U(r) = \int_{r_0}^r p(v(r)) \frac{dv}{dr} dr = \int_{r_0}^r p(\pi
  r^2) 2\pi r dr.
\end{equation*}
That work is used to deform and accelerate the cladding.  We say that
the energy of deformation is
\begin{equation*}
  D(r) = \sigma \log \left(\frac{r}{r_0} \right),
\end{equation*}
and the kinetic energy of the cladding is
\begin{equation*}
  K = \frac{1}{2} m \left( \dot r \right)^2.
\end{equation*}
Thus the velocity as a function of radius is
\begin{equation*}
  v(r) = \sqrt{\frac{2}{m} \left( U(r) - \sigma \log
      \left(\frac{r}{r_0} \right) \right) },
\end{equation*}
and we can obtain the velocity as a function of time by numerical
integration.

\section{Splines}
\label{sec:splines}

From Fig.~\ref{fig:basis}, I conclude that the second derivative is
affine between knots and that consequently if it is positive at the
knots it is positive between them.  Further, if the second derivative
is positive over the domain and the first derivative is positive at
the right hand end point, the the function is monotonic and convex
over the entire domain.
\begin{figure}
  \centering
    \resizebox{\columnwidth}{!}{\includegraphics{basis.pdf}}  
    \caption{Cubic spline basis functions and their first and second
      derivatives. The dashed plots are for $f_5$ whose knot is
      $t_5=3$.  Note that $t_k=0$ for $f_0,~f_1,~f_2$ and $f_3$, and
      $t_4=2$.}
  \label{fig:basis}
\end{figure}

Figure~\ref{fig:basis} also illustrates the symmetry of the right and
left boundaries.  The knot location $t$ of each basis function is the
point on the left where the function departs from zero.  There are
four functions at $t=0$.  There are also four functions at $t=10$, but
since their coefficients are always zero, they are not degrees of
freedom in any fit.

\section*{Appendix}
\label{sec:appendix}

See other information in cmfSuite/doc/HE.tex

One may fetch source code for this project from Google code and build
this document as follows:
\begin{verbatim}
->git clone https://code.google.com/p/metfie
->cd metfie/like_gun
->make notes.pdf
\end{verbatim}

%
\vfill \hrule

Source file: https://code.google.com/p/metfie/like\_gun/notes.tex

\end{document}

%%%---------------
%%% Local Variables:
%%% eval: (TeX-PDF-mode)
%%% eval: (setq ispell-personal-dictionary "./localdict")
%%% End:
